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1  Introduction 


Background 


HrVEL2D  is  a  free-surface,  depth-averaged,  two-dimensional  model 
designed  specifically  for  flow  fields  that  contain  supercritical  and  subcritical 
regimes  as  well  as  the  transitions  between  the  regimes.  The  model  provides 
numerically  stable  solutions  of  advection-dominated  flow  fields  containing 
shocks  such  as  oblique  standing  waves  and  hydraulic  jumps. 

HrVEL2D  has  been  verified  by  comparing  computed  model  results  with 
laboratory  data.  The  findings  of  these  tests  are  presented  in  Stockstill  and 
Berger  (1994). 


Purpose  and  Scope 


This  report  serves  as  a  user’s  manual  for  HIVEL2D  Version  1.07.  First,  a 
brief  description  of  the  model  equations  are  presented.  Next,  step-by-step 
instructions  are  given  to  guide  the  user  in  creating  the  required  input  files. 
These  input  files  include  the  geometry  (numerical  model  computational  mesh) 
file,  the  hydrodynamic  parameters/boundary  conditions  file,  and  the  initial 
conditions  file.  Finally,  the  numerical  model  output  file  is  explained. 
Appendix  A  describes  the  boundary  conditions.  Appendix  B  describes 
problems  with  the  model  and  solutions.  Appendix  C  gives  an  example 
problem.  Appendix  D  discusses  the  interpolation  program. 
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Introduction 


2  HIVEL2D  Overview 


HIVEL2D  (Stockstill  and  Berger  1994)  is  designed  to  simulate  flow  typical 
in  high-velocity  channels.  The  model  is  a  finite  element  description  of  the 
two-dimensional  shallow- water  equations  in  conservative  form.  The  model 
does  not  include  Coriolis  or  wind  effects  as  these  are  typically  not  important 
in  high-velocity  channels. 


Governing  Equations 

Vertical  integration  of  the  equations  of  mass  and  momentum  conservation 
for  incompressible  flow  with  the  assumption  that  vertical  accelerations  are 
negligible  compared  to  the  acceleration  of  gravity  results  in  the  governing 
equations  commonly  referred  to  as  the  shallow-water  equations.  The 
dependent  variables  of  the  two-dimensional  fluid  motion  are  defined  by  the 
flow  depth  h,  the  volumetric  discharge  per  unit  width  in  the  x-direction  p,  and 
the  volumetric  discharge  per  unit  width  in  the  y-direction  q.  These  variables 
are  functions  of  the  independent  variables  x  and  y,  the  two  space  directions, 
and  time  t.  The  shallow-water  equations  in  conservative  form  are  given  as 
(Abbott  1979): 


dt  dx  dy 


+  =  0 


(1) 


where 


Q  = 


\  ^  ) 


(2) 
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where 

p  =  uh,  u  being  the  depth-averaged  x-direction  component  of 
velocity 

q  =  vh,  V  being  the  depth-averaged  y-direction  component  of 
velocity 

g  =  acceleration  due  to  gravity 

a^,ay^,ay^,ayy  =  Reynolds  stresses  where  the  first  subscript  indicates  the 

direction  and  the  second  indicates  the  axis  direction  normal 
to  the  face  on  which  the  stress  acts 


p  =  fluid  density 


Zg  =  chaimel  invert  elevation 


n  =  Manning’s  roughness  coefficient 
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Cq  =  dimensional  constant  {€„  =  1  for  SI  units  and  Q=  1.486 
for  non-SI  units) 

The  Reynolds  stresses  are  determined  using  the  Boussinesq  approach  relating 
stress  to  the  gradient  in  the  mean  currents: 


=  2pv , 


dll 


a  =  o  =  pv, 

Jey  yr  r  I 


dx 
(  du 


dv\ 


dy  dx 


o  =  2pv, _ 

^  'ay 


(6) 


where  v,  is  the  eddy  viscosity,  which  varies  spatially  and  is  solved  empirically 
as  a  function  of  local  flow  variables  (Rodi  1980); 


V,  = 


Cn 


^ ii 


(7) 


where  C  is  a  coefficient  that  varies  between  0.1  and  1.0. 

The  typical  high-velocity  channel  is  rather  wide  and  shallow,  so  that  the 
bulk  of  the  turbulence  is  bed-generated  and  therefore  Equation  7  is  reasonable. 
However,  if  the  walls  or  other  structures  make  a  significant  contribution  to  the 
overall  turbulence,  this  equation  will  underpredict  the  magnitude  of  the  eddy 
viscosity  or  turbulence.  Even  under  typical  channel  configurations  there  is 
fairly  large  uncertainty  in  the  parameter  C.  This  is  apparent  since  it  has  been 
found  to  range  from  0.1  to  1.0.  Therefore,  the  user  should  consider  making  a 
sensitivity  check  on  the  calculated  water  surface  for  a  range  of  values  of  C  and 
any  other  parameters  having  values  that  are  largely  uncertain. 


Finite  Element  Model 

This  system  of  partial  differential  equations  is  solved  using  the  finite  ele¬ 
ment  method.  The  finite  element  approach  taken  is  a  Petrov-Galerkin  formula¬ 
tion  that  incorporates  a  combination  of  the  Galerkin  test  function  and  a 
non-Galerkin  component  to  control  oscillations  due  to  convection.  An  integra¬ 
tion  by  parts  procedure  is  used  to  develop  the  weak  form  of  the  equations. 

The  weak  form  facilitates  the  specification  of  boundary  conditions.  It  is  given 
as: 
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V  ‘  dt  ax  '  dy  ^  ‘  dx 

+  (p  +'^H  dO.^  +  /4>i  ■*■  ~  ® 

'  dy  '  ‘  i 


where  the  variables  are  understood  to  be  discrete  values  and 


e  =  subscript  indicating  a  particular  element 


Q  =  domain 


ip.  =  (J)/  +  <Pi  =  test  function 
(f);  =  Galerkin  part  of  the  test  function 
I  =  identity  matrix 

cp;  =  non-Galerkin  part  of  the  test  function 
(n^  n  )  =  n  =  unit  vector  outward  normal  to  the  boundary  F, 


Natural  boundary  conditions  are  applied  to  the  sidewall  boundaries  through  the 
weak  form.  The  sidewall  boundaries  are  “no-flux”  boundaries;  that  is,  there  is 
no  net  flux  of  mass  or  momentum  through  these  boundaries.  These  boundary 
conditions  are  enforced  through  the  line  integral  in  the  weak  form. 


Petrov-Galerkin  Test  Function 


The  Petrov-Galerkin  test  function  ip,  is  defined  (Berger  1993)  as: 


'H’i  =  4>/  +  «Pi 
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where 


9;  =  P 


/  \ 
3(j) .  «  ^(f) .  - 

Ax—IA  +  Ay  ‘B 
dx  dy  j 


(11) 


where  P  is  a  dissipation  coefficient  varying  in  value  from  0  to  0.5,  the  4)  terms 
are  the  linear  basis  functions,  and  Ajc  and  Ay  are  the  ^rid  intervals.  A  detailed 
explanation  of  this  test  function,  in  particular  A  and  B,  is  given  in  Berger 
(1993). 

Shock  Capturing 

The  coefficient  P  scales  the  dissipation  needed  for  numerical  stability. 

More  dissipation  is  needed  in  the  vicinity  of  shocks  such  as  hydraulic  jumps 
than  in  smooth  regions  of  the  flow  field.  Because  a  lower  value  of  p  (P  = 
is  more  precise,  a  large  value  of  p  (P  =  Pjy,  =  0.5  where  P^;^  and  P^^  are 
the  Petrov-Galerkin  parameters  for  smooth  flow  and  for  shocks,  respectively) 
is  applied  only  in  regions  in  which  it  is  needed.  HIVEL2D  employs  a  mecha¬ 
nism  that  detects  shocks  and  increases  P  automatically.  Therefore,  P^^/  is 
implemented  only  when  needed  as  determined  by  evaluation  of  the  element 
energy  deviation.  In  a  similar  manner,  the  eddy  viscosity  coefficient  C  varies 
from  CsM  to  C^,  the  effect  being  that  eddy  viscosity  is  increased  only  in  areas 
of  greatest  element  energy  deviation. 


Temporal  Derivatives 

A  finite  difference  expression  is  used  for  the  temporal  derivatives.  The 
general  expression  for  the  temporal  derivative  of  a  variable  Qj  is: 


/  ^ 
dQj 

m+1 

'  1  +  a' 

'qT'  -  Qf 

\  J 

+ 

i  2  J 

fl  -  a] 

/m  +  1  f  m 

'q:  -  QT'' 

K  2  J 

^  _  f-i 

(12) 


where  j  is  the  nodal  location  and  m  is  the  time-step.  An  a  equal  to  1  results 
in  a  first-order  backward  difference  approximation,  and  an  a  equal  to  2  results 
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in  a  second-order  backward  difference  approximation  of  the  temporal 
derivative. 


Solution  of  the  Nonlinear  Equations 

The  system  of  nonlinear  equations  is  solved  using  the  Newton-Raphson 
iterative  method.  Let  J?,  be  a  vector  of  the  nonlinear  equations  computed  using 
a  particular  test  function  and  using  an  assumed  value  of  Qj.  is  the  resid¬ 
ual  error  for  a  particular  test  function  i.  Subsequently,  /?,  is  forced  toward  zero 
as: 


(13) 


where  the  derivatives  composing  the  Jacobian  are  determined  analytically  and 
k  is  the  iteration  number.  This  system  of  equations  is  solved  for  AQj  and  then 
an  improved  estimate  for  Qf*’  is  obtained  from: 


=  Qj  +  Afi* 


(14) 


This  procedure  is  continued  until  convergence  to  an  acceptable  residual  error  is 
obtained. 


Model  Features 

Particular  features  of  HIVEL2D  Version  1.07  are  as  follows: 

a.  Combinations  of  linear-based  triangular  and  rectilinear  shape  functions 
are  used  to  represent  p,  q,  and  h. 

b.  The  model  is  compiled  to  solve  problems  having  meshes  as  large  as 
2000  nodes  and  2000  elements  and  will  run  on  a  386-based  (or  higher) 
IBM-compatible  personal  computer  having  8  MB  of  RAM. 

c.  A  Newton-Raphson  iteration  is  used  within  each  time-step. 

d.  The  mild  slope  assumption  is  invoked.  This  means  the  slope  should  be 
less  than  about  0.05.  This  assumes  that  the  slope  is  geometrically  mild, 
not  hydraulically  mild. 
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e.  Boundary  conditions  are  specified  for  any  combination  of  supercritical 
and  subcritical  inflow  and  outflow. 

/.  Boundary  conditions  are  constant  over  the  simulation  period. 

g.  Boundary  roughness  (Manning’s  n)  is  specified  on  an  element  type 
basis. 

h.  A  user-specified  parameter  a  produces  either  first-order  or  second-order 
backward  temporal  differences. 

i.  A  Petrov-Galerkin  approach  in  which  the  test  function  is  weighted 
upwind  along  characteristics  is  employed.  The  degree  of  upwinding 
and  thus  stability  is  determined  by  the  parameters  and  jS^^. 

Defaults  are  =  0.25  and  =  0.5. 

j.  A  shock  detection  mechanism  based  upon  the  energy  variation  per 
element  is  used  to  invoke 

k.  Turbulent  eddy  viscosity  is  calculated  based  upon  simple  user-specified 
parameters  and  using  velocity,  depth,  and  roughness. 
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Developing  an  Application 


The  basic  steps  to  developing  an  application  of  HIVEL2D  are  as  follows: 

a.  Generate  the  grid: 

(1)  Number  the  grid  intelligently. 

(2)  Identify  inflows  and  outflows. 

b.  Develop  a  hydrodynamic  input  file. 

c.  Develop  a  reasonable  hot  start  file. 

d.  Run  the  model  (probably  several  times).  If  necessary,  refine  the  grid 
and  interpolate  a  new  hot  start  file  using  INTRPL8.EXE  (Appendix  D). 

e.  Run  HIVELBIN.EXE  on  the  output  file  to  generate  a  solution  file  that 
can  be  viewed  in  FastTABS. 

f.  Examine  the  solution  for  reasonableness. 

Grid  Generation 

To  a  large  degree  the  quality  of  the  grid  determines  the  accuracy  and 
stability  of  the  model.  The  first  step  in  grid  generation  is  getting  the  necessary 
geometry  information  into  the  grid  generator.  Critical  elements  such  as  points 
on  transitions,  bridge  piers,  and  curves  should  be  put  into  the  grid  generator. 
One  method  involves  drawing  the  structure  to  be  modeled  in  three  dimensions 
in  a  CAD  program,  and  then  extracting  the  critical  points.  These  critical 
points  are  then  put  into  a  grid  generator.  One  grid  generator  the  Corps  of 
Engineers  uses  is  FastTABS  (Brigham  Young  University  1993).  FastTABS  is 
a  pre-  and  postprocessor  for  programs  involving  two-  dimensional  finite  ele¬ 
ment  meshes.  It  was  developed  by  the  Engineering  Computer  Graphics 
Laboratory  at  Brigham  Young  University  for  the  U.S.  Army  Engineer 
Waterways  Experiment  Station. 
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Remember  the  following  general  rules  when  generating  the  grid: 


a.  HIVEL2D  uses  linear  elements.  Therefore,  when  generating  a 
HIVEL2D  grid,  use  only  four-node  quadrilaterals  and  three-node 
triangles. 

b.  Keep  the  element  aspect  ratio  less  than  3:1  (the  closer  to  1:1,  the  better). 
The  aspect  ratio  is  the  ratio  of  the  longest  element  dimension  to  the 
shortest,  i.e.,  the  length-to-width  ratio. 

c.  Use  gradual  transitions  in  element  size.  Generally,  an  element’s  area 
should  not  be  greater  than  IVi  times  its  smallest  neighbor. 

d.  Include  at  least  five  or  six  elements  across  a  channel  in  the  area  of 
interest.  If,  for  example,  the  channel  has  an  island  in  the  center,  this 
resolution  is  needed  on  both  sides  of  the  island.  Also,  increase  the 
resolution  around  grade  breaks  and  wall  transitions. 

Once  the  grid  generation  is  complete,  be  sure  to  renumber  the  grid  using 
FastTABS’  renumber  options.  The  best  numbering  scheme  will  give  the 
smallest  bandwidth.  Normally,  the  nodes  need  to  be  numbered  progressively 
across  the  narrowest  dimension  of  the  grid.  This  minimizes  the  bandwidth, 
which  makes  HIVEL2D  run  more  efficiently. 

Before  exiting  the  grid  generator,  note  the  nodes  that  compose  the  inflow 
and  outflow  boundaries.  These  nodes  must  be  put  into  the  geometry  (grid) 
file.  Open  the  geometry  file  with  an  ASCII  text  editor,  such  as  the  MS-DOS 
editor  EDIT.COM.  Add  the  BI  lines  for  each  node  in  which  an  inflow  condi¬ 
tion  is  to  be  input.  Add  a  BO  line  for  each  outflow  node  (this  will  not  always 
require  a  boundary  condition).  These  lines  are  no  longer  formatted  in  v.1.07 
of  HIVEL2D.  Be  sure  to  type  BO,  not  BO  (zero).  Keep  in  mind  that  the 
model  is  expecting  a  relatively  narrow  (computationally)  dimension  across  the 
channel.  The  various  lines  can  appear  in  any  order.  Any  prefixes  on  each  line 
that  are  not  required  by  HIVEL2D  are  printed  to  the  screen,  but  are  otherwise 
ignored  by  the  program. 

Table  1  shows  the  geometry  file  format.  If  the  geometry  was  generated 
using  FastTABS,  then  the  GNN  and  GE  lines  are  transparent  to  the  user.  To 
generate  a  grid  from  scratch,  be  sure  to  follow  the  correct  format  shown  in 
Table  1. 


Hydrodynamic  Input 

Starting  the  model  is  usually  the  most  difficult  task.  The  model  should  be 
started  with  small  time-steps  (i.e.,  a  CFL  number  less  than  1)  where: 
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Table  1 

Geometry  File  Format 

Prefix 

(Columns  1-3) 

Column 

Read  As.. 

Contents 

Bl  or  bi 

Free  format 

Node  on  the  boundary  in  which  flow  enters  the 
model 

BO  or  bo 

Free  format 

Node  on  the  boundary  in  which  flow  leaves  the 
model 

GNN  or  gnn^ 

Free  format 

Node  number, 

node's  x-coordinate, 

node's  y-coordinate, 

node's  z-coordinate  (invert  elevation) 

GE  or  ge^ 

4-8 

Element  number 

9-14 

Node  number  1  of  element  number  GE 

15-20 

Node  number  2  of  element  number  GE 

21-26 

Node  number  3  of  element  number  GE 

27-32 

Node  number  4  of  element  number  GE 
(0  if  element  is  a  triangle) 

56-61 

Material  tvoe  -  used  to  designate  the  bed  roughness 
area  of  common  type.  The  actual  value  of 

Manning's  n  is  input  in  the  hydrodynamic  input  file. 

The  user  can  specify  maximum  of  three  material 
types. 

’Need  one  GNN  line  for  each  node.  The  GNN's  do  not  have  to  be  in  order,  but  no  node 
numbers  can  be  left  out  between  1  and  the  largest  node  number  used. 

^Need  one  GE  line  for  each  element  giving  the  list  of  nodes  associated  with  this  element 
in  counterclockwise  order. 

CFL  = 

Al  /At 

(15) 

where 

I  U\  =  Vm^+ 


and 


A/  =  the  element  length 
At  =  the  time-step  size 
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The  time-step  can  then  be  gradually  increased.  If  the  steady-state  solution  is 
desired,  then  a  fairly  large  time-step  can  be  used.  If  the  interest  is  in 
unsteady  results,  then,  for  accuracy,  the  CFL  number  should  be  held  to  a 
maximum  of  2. 

Appropriate  boundary  conditions  must  be  supplied.  For  an  inflow  bound¬ 
ary  condition  when  the  flow  is  supercritical,  the  user  must  specify  both  x- 
and  y-components  of  flow  along  with  the  depth.  If  the  flow  is  subcritical 
here,  then  only  the  x-  and  y-components  of  flow  are  read.  If  the  flow  is  sub- 
critical  at  the  outflow  boundary,  then  the  tailwater  elevation  must  be 
specified.  If  the  flow  is  supercritical  at  the  outflow,  then  no  boundary  con¬ 
ditions  need  to  be  given.  This  is  done  in  the  input  file  by  specifying  a 
tailwater  elevation  that  is  lower  than  the  bed  elevation  at  the  outflow 
boundary. 

Starting  the  model  is  sometimes  difficult  when  the  flow  is  supposed  to  be 
supercritical  at  the  outflow  boundary.  The  model  attempts  to  converge  to  a 
solution,  but  it  may  not  be  the  desired  solution.  One  way  to  avoid  this  is  to 
start  the  model  with  a  tailwater  that  is  slightly  subcritical  (if  the  starting  con¬ 
ditions  are  not  known).  Then  after  the  model  settles  down,  the  boundary 
condition  can  be  changed  to  supercritical.  This  takes  a  little  experience. 

Users  should  note  that  the  Manning’s  n  applies  to  each  element  type  as 
well  as  the  adjoining  sidewalls.  It  is  also  used  in  conjunction  with  velocity, 
depth,  and  C  to  determine  a  turbulent  eddy  viscosity  estimate. 

The  input  to  the  hydrodynamics  file  is  made  through  the  list  in  Table  2. 
All  variables  are  real  except  those  beginning  with  /,  J,  K,  L,  M,  or  N,  which 
are  integer  variables.  All  input  is  free-field.  Each  variable  is  described  in 
Table  2  in  the  order  in  which  it  appears  in  the  hydrodynamic  input  file. 

Units  are  designated  by  L  and  T,  for  length  and  time,  respectively.  If  no 
units  are  specified,  the  variable  is  dimensionless. 


Hot  Start  File 

This  file  is  the  last  two  time-steps  of  information  from  the  previous  model 
run.  Since  the  temporal  derivative  can  be  a  second-order  backward  differ¬ 
ence,  two  time-steps  of  old  information  are  needed.  If  this  is  the  initial  run, 
then  this  information  will  have  to  be  generated  based  upon  experience.  One 
method  is  to  start  the  model  with  zero  velocity  and  a  constant  depth.  The 
data  are  read  as  free-field. 

This  file  will  be  overwritten  by  the  model  run,  so  if  there  is  a  chance  that 
this  information  will  be  needed  later,  save  a  copy.  If  the  model  run  finishes 
but  the  data  are  bad  for  some  reason,  the  hot  start  will  be  bad  as  well. 

Again,  this  is  a  good  reason  to  save  a  copy. 
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Table  2 

Hydrodynamic  Input  File  Format 

Line  Number 

Variables 

Description 

1 

9 

Gravitational  acceleration  (REAL),  i/T* 

2 

C\  ( C2 ) 

(REAL)  is  the  empirical  conversion  coefficient  used 
in  Manning's  equation.  This  needs  to  be  consistent 
with  the  units  used  for  g.  Typical  values  are  2.208  for 
English  units  and  1 .0  for  SI. 

3 

These  coefficients  (REAL)  supply  the  turbulence 

(CSB.CSBSHK) 

coefficients  in  the  normal  flow  region  and  near  the 
jump,  respectively.  0.1s  Cj^yS  1.0  Typically, 

is  chosen  as  0.1  and  Cg„  is  0.5.  These  parameters 
determine  the  turbulent  eddy  viscosity  based  upon  the 
depth,  velocity,  and  roughness.  Larger  values  imply 
larger  turbulent  eddy  viscosity.  Again,  if  the  flow  is 
relatively  smooth  and  there  are  no  jumps,  set  Cg^  to 
the  value  used  for  Cg^ 

4 

At,  a 

At  (REAL)  is  the  time-step  size  T  (units  need  to  be 

(DELT, 

consistent  with  other  variables). 

ALPHAT) 

a  (REAL)  determines  the  precision  of  the  temporal  term 
in  the  model. 

1.0  s  a  s  2.0 

a  =  1 .0  for  first-order  backward  difference 
a  =  2.0  for  second-order  backward  difference 

On  steady-state  problems  use  1.0  fit  is  more  stable). 

On  time-dependent  problems  generally  use  2.0,  but  if 
the  time-step  is  small  it  does  not  matter. 

5 

NSTEPS,  NRVL 

NSTEPS  (INTEGER)  is  the  number  of  time-steps  to 
model  in  this  run.  NRVL  (INTEGER)  is  the  interval  for 
which  all  hydrodynamic  information  is  written  to  the  hy¬ 
drodynamic  output  file,  which  is  for  plotting.  Small  val¬ 
ues  for  this  can  cause  the  output  file  to  become  quite 
large.  If  interest  is  only  on  the  last  time-step,  set 
NRVL=NSTEPS. 

For  example,  entering  a  value  of  2  will  save  fme-steps 

2,  4.  6.... 

(Sheet  1  of  3) 
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Table  2  (Continued) 

Line  Number 

Variables 

Description 

6 

UMAX,  TOL 

UMAX  (INTEGER)  is  the  maximum  number  of  Newton- 
Raphson  iterations  per  time-step.  A  typical  value  is  4. 
TOL  (REAL)  is  the  convergence  criteria.  The  number 
of  Newton-Raphson  iterations  per  time-step  is  UMAX 
unless  convergence  is  reached  first.  A  typical  value  is 
0.005. 

If  ]  1  AA  1 1  <  TOL,  then  the  solution  is  considered 
converged  and  no  more  iterations  are  performed  for 
the  current  time-step: 

Where 

1 1  AA  1 1  =  maximum  of  |  (  A  -  Aq  )  ]  /  A  (for  each 
node) 

A  =  Froude  number  this  iteration 

Fq  =  Froude  number  last  iteration 

7 

lUORPj.ND,, 

VXi,  VY;,  Hi 

These  are  the  inflow  input  lines.  One  line  is  supplied 
for  each  inflow  node  (these  are  the  nodes  identified  by 

Bl  in  the  geometry  file). 

ND|  (INTEGER)  is  the  individual  node  number. 

VX|,  VYj  (REALS)  are  the  components  of  either  veloc¬ 
ity  (w,  v),  L/T,  or  unit  discharge  (p,  q),  L^fT,  at  this 
node.  The  type  is  determined  by  IUORP|  (Table  3). 

H,  (REAL),  L,  is  the  node  depth.  This  is  used  for 
supercritical  inflow  boundary  conditions.  H  j  is  used  or 
ignored  depending  on  lUORPi  (Table  3).  It  is  not 
necessary  to  put  a  value  here  if  W,-  is  to  be  ignored. 

8 

TAIL 

TAIL  (REAL)  is  the  tailwater  elevation,  L,  (this  will  be 
assigned  to  all  nodes  that  have  been  designated  as 
outflow  nodes,  BO  lines,  in  the  geometry  file).  If  this 
elevation  is  less  than  the  bed  elevation  at  this  node, 
the  model  assumes  that  the  flow  is  supercritical  there 
and  no  boundary  condition  is  assigned. 

9 

NMTYP 

NMTYP  (INTEGER)  is  the  number  of  element  types  for 
roughness.  This  list  begins  with  type  1  and  does  not 
skip  an  integer.  These  element  types  are  assigned  in 
the  geometry  file. 

(Sheet  2  of  3) 
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Table  2  (Concluded) 


Line  Number 

Variables 

Description 

10 

j,  FRICj 

j  (INTEGER)  is  the  element  type. 

FRIC  J  (REAL)  is  the  Manning’s  n  roughness  for 
element  type  j.  There  should  be  NMTYP  number  of 
these  lines,  beginning  with  j=1  through  j=NMTYP,  in 
order,  with  no  skips.  This  friction  applies  to  the  bed 
and  the  sidewalls  of  the  elements  of  this  type. 

Note:  Ttie  parameters  and  (Petrov-Galerkin  weight  coefficients)  are  not  required  in 
the  hydrodynamic  input  file.  The  defaults  are  -  0.25  and  =  0.5.  These  values  can  be 
overridden  by  including  as  the  first  line  of  the  hydrodynamic  input  (for  example): 


OR  0.20  0.5 

where  OR  is  the  flag  that  Instructs  the  program  to  read  the  following  numbers 
as  =  0.20  and  /Sg^  =  0.5 


(Sheet  3  of  3) 


Table  3 

Designation  of  Boundary  Condition  Type 

Flow  Regime 

Specified  Variables 

lUORP 

Supercritical 

P,  q.  h 

2 

Supercritical 

u,  V,  h 

-2 

Subcritical 

p.  q 

1 

Subcritical 

U,  V 

-1 

The  information  in  the  hotstart  file  is: 

(Line  1)  TIME 

(Line  2)  q^,  ff,  cT'^ ,  (for  each  node,  in  order) 


where 

TIME  =  the  time  associated  with  values  at  time-step  m 
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m  =  the  last  time-step 
m-1  =  the  next  to  the  last  time-step 
Remember  that  this  file  is  free  format. 

The  initial  hot  start  file  can  also  cause  starting  a  run  to  be  tricky.  If  the 
interest  is  in  steady-state  conditions,  the  accuracy  of  the  initial  guess  will 
determine  how  long  it  wilt  take  to  reach  steady  state.  If  the  geometry  is 
simple,  then  there  is  no  trouble  in  guessing  the  discharge  components.  How¬ 
ever,  if  it  is  complicated,  it  may  be  necessary  to  resort  to  zero  discharge  and  a 
specified  depth,  which  is  a  typical  method.  If  the  walls  are  sloping,  there 
may  be  some  difficulty  in  specifying  the  depth  as  well.  It  is  fairly  simple  to 
write  a  program  to  subtract  the  bed  elevation  from  some  set  of  water-surface 
elevations.  Typically,  the  user  should  set  the  old  and  the  even  older  time- 
step  information  in  this  file  as  identical. 


Running  the  Model 

To  run  the  model  the  user  must  supply  four  file  names,  three  input  and 
one  output.  Note  the  following  example  (the  computer  prompts  are  in  bold 
and  the  user’s  commands  and  responses  are  in  italics): 

C-\>HIVEL2D 

WHAT  IS  THE  GEOMETRY  FILE  ? 
example.geo 

WHAT  IS  THE  HYDRODYNAMIC  AND  COMPUTATION 
PARAMETER  INPUT  FILE  ? 

example.flo 

WHAT  IS  THE  PLOT  OUTPUT  FILE  ? 
example. out 

WHAT  IS  THE  HYDRODYNAMIC  OUTPUT  FILE  (THIS  IS 
ALSO  THE  NAME  OF  THE  INITIAL  DATA  FILE  FOR 
HOTSTART)  ? 
example.hot 

a.  The  executable  is  called  HIVEL2D.EXE. 

b.  Each  file  name  can  be  up  to  15  characters  for  Unix  and  Apple  machines, 
and  12  characters  for  IBM  compatibles  (e.g.,  EXAMPLEl.GEO). 

c.  The  file  extensions  of  GEO,  FLO,  OUT,  and  HOT  are  ones  that  are 
typically  used  to  designate  files,  but  these  suffixes  are  arbitrary. 

d.  The  geometry  file  contains  the  node  locations,  the  element  connection 
table,  and  the  designation  of  boundary  inflow  and  outflow  nodes. 
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e.  The  hydrodynamic  and  computation  parameter  input  file  contains  the 
hydraulic  information  about  the  run  as  well  as  the  computation 
parameters.  Examples  of  computation  parameters  are  the  specific 
boundary  values  of  depth  and  flows,  roughness,  time-step  size,  number 
of  iterations,  etc. 

/.  The  output  file  contains  the  water-surface  elevation  and  velocity  fields 
for  all  nodes  at  each  time  location  designated  in  the  hydrodynamic  input 
file.  This  information  can  then  be  converted  to  a  form  that  can  be 
viewed  in  a  postprocessor  such  as  FastTABS. 

g.  The  hot  start  file  contains  the  time,  p,  q,  and  depth  for  the  last  two  time- 
steps.  This  is  used  to  start  the  calculation.  This  file  is  overwritten  by 
each  successful  run;  so,  if  it  might  be  needed  later,  make  a  copy  of  it 
under  another  name  before  running  the  program. 

After  the  file  names  are  input,  the  program  will  begin  running  and  several 
banners  will  appear  on  the  screen. 

a.  First,  the  program  informs  the  user  of  any  unused  lines  in  the  geometry 
file  and  prints  the  contents  of  the  unused  lines. 

b.  Next,  the  contents  of  the  hydrodynamic  input  file  are  displayed  on  the 
screen  to  allow  the  user  to  check  the  values  that  were  read. 

c.  From  the  initial  conditions  the  average  energy,  the  minimum  vorticity 
and  the  element  where  it  occurred,  and  the  maximum  vorticity  and  the 
element  where  it  occurred  are  calculated  and  displayed. 

d.  If  the  coefficients  C^,  C^h,  or  a  are  out  of  range,  the  program  will  re¬ 
port  the  error  and  halt. 

e.  Then,  as  the  program  runs,  the  results  for  each  time-step  are  displayed. 
These  results  include  the  number  of  iterations  required,  the  maximum 
residual  error  and  the  node  with  which  it  is  associated,  the  average  ener¬ 
gy,  the  minimum  vorticity  and  the  element  where  it  occurred,  and  the 
maximum  vorticity  and  the  element  where  it  occurred. 


Output  File 

This  file  is  output  from  the  HIVEL2D  run  and  is  intended  mainly  to  supply 
information  for  plotting  in  FastTABS.  The  output  file  contains  the  following: 

a.  The  number  of  time-steps  saved,  number  of  elements,  and  number  of 
nodes. 

b.  The  time. 
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c.  Velocity  components  of  u  and  v,  and  water-surface  elevation  for  each 
node. 

d.  Material  type  number  for  each  element. 

Note:  these  are  repeated  for  each  time-step  that  was  output. 

The  number  of  time-steps  saved  is  determined  by  the  request  in  the  hydro- 
dynamic  input  file  for  calculated  time-steps  and  NRVL  (the  interval  to  save 
time-steps  to  this  file).  This  file  can  get  quite  large,  so  assign  NRVL  with 
discretion. 


Viewing  the  Output 

The  output  files  generated  by  HIVEL2D  are  ASCII  text  files,  which,  with 
some  manipulation,  can  be  imported  to  almost  any  numerical  postprocessor 
available.  The  preferred  postprocessor  is  FastTABS,  because  the  output  files 
are  in  a  format  that  can  be  directly  converted  to  the  appropriate  binary  form 
required  by  FastTABS.  The  program  that  converts  HIVEL2D  output  files  to 
FastTABS  binary  solution  files  is  called  HIVELBIN.  The  required  version  of 
HIVELBIN  depends  on  the  FastTABS  version  used. 
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Appendix  A 
Boundary  Conditions 


Flow  In  and  Out 

Case  1 


Flow  regime 
Supercritical  upstream 
Supercritical 
downstream 

Boundary  conditions 
p,  q,  h  specified 
upstream  (defined  in 
Chapter  2) 

Comments 

It  can  be  difficult  to  get 
the  outflow  boundary  to  converge  at  startup.  One  method  is  to  specify  a 
tail  water  such  that  the  flow  is  barely  subcritical,  then  after  some  time,  when 
the  resulting  jump  is  nearly  gone,  remove  the  downstream  boundary 
condition.  Note:  F,  =  inflow  Froude  number;  F^  =  outflow  Froude  number. 

Case  2 


Flow  regime 
Supercritical  upstream 
Subcritical  downstream 

Boundary  conditions 

p,  q,  h  specified 
upstream 

h  specified  downstream 

Figure  A2.  Case  2  sketch 
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Case  3 


Flow  regime 
Subcritical  upstream 
Supercritical  downstream 

Boundary  conditions 
p,  q  specified  upstream 


Figure  A3.  Case  3  sketch 

Case  4 

Flow  regime 
Subcritical  upstream 
Subcritical  downstream 

Boundary  conditions 

p.  q  specified  upstream 
h  specified  downstream 

Comments 

This  problem  may  take  a  Figure  A4.  Case  4  sketch 
long  time  to  reach 
steady  state. 


Case  5 

Flow  regime 
Subcritical  upstream 
Subcritical  downstream 

Boundary  conditions 

p,  q  specified  upstream 
h  specified  downstream 

Comments 
Even  though  this  con¬ 
tains  a  steep  section  in 
which  the  flow  is  supercritical,  the  boundary  conditions  are  the  same  as  for 
Case  4.  The  tailwater  sfiecified  downstream  causes  the  hydraulic  jump. 


Figure  A5.  Case  5  sketch 
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Flow  In,  No  Flow  Out 


Case  6 

Flow  regime 
Supercritical  upstream 
Wall  downstream 

Boundary  conditions 
p,  q,  h  upstream 
Leave  out  BO  lines 
(BI  and  BO  defined  in 
Table  1) 


Figure  A6.  Case  6  sketch 


Comments 

No  steady  state,  keeps  filling.  Note  the  moving  jump  that  forms. 


Case  7 

Flow  regime 
Subcritical  upstream 
Wall  downstream 

Boundary  Conditions 

p,  q  specified  upstream 
Leave  out  BO  lines  in 
geo  file 


Comments  Figure  A7.  Case  7  sketch 

There  is  no  steady- 
state  answer,  just 
keeps  filling. 


No  Flow  In,  Flow  Out 


Case  8 

Flow  regime 
Wall  upstream 
Supercritical  downstream 

Boundary  conditions 
Leave  out  BI  lines 
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Comments 


No  steady  state,  goes  dry.  This  is  probably  a  very  hard  problem  to  run. 


Case  9 

Flow  regime 

Wall  upstream 
Subcritical  downstream 

Boundary  conditions 
Leave  out  Bl  lines 
Specify  h  downstream 

Comments 

No  steady  state,  goes 
dry.  This  is  also  a  very 
hard  problem  to  run. 


No  Flow  In  or  Out 

Case  10 


Flow  regime 

Wall  upstream 
Wall  downstream 

Boundary  conditions 
Leave  out  Bl  and  BO 
lines 


Figure  A10.  Case  10  sketch 
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Appendix  B 
Troubleshooting 


Table  B1 

Problems  and  Solutions 

Problem 

Possible  Solutions 

1 .  Model  won't  converge 

Try  more  iterations  per  time-step. 

Reduce  time-step. 

If  lower  boundary  is  migrating  away  from 
solution;  First  set  a  tailwater,  wait  until  jump 
is  nearly  gone,  then  remove  tailwater 
boundary  condition. 

Call  for  help;  be  prepared  to  send  the 
geometry,  hydraulic  input,  output,  and  hot 
start  files. 

2.  Model  keeps  filling  up  (water  level 
continuously  rises). 

Make  sure  you  didn't  leave  out  BO  lines  in 
geometry  file.  Use  BO  and  not  BO  (zero) 

(Table  1). 

3.  Converges  to  a  bizarre  solution. 

Make  sure  resolution  is  adequate. 

4.  Problems  reading  hot  start  file. 

Verify  that  the  time  is  included  as  the  first 
line  of  the  file. 

Verify  the  number  of  lines. 

May  have  to  put  a  carriage  return  at  the  end 
of  the  last  line  of  the  hot  start  file. 
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Appendix  C 
Example  Problem 


The  example  grid  is  shown  in  Figure  Cl.  The  flow  enters  on  the  left  and 
exits  to  the  right.  The  2,000-ft  (610-m)-long  model  is  supercritical  at  both 
boundaries,  so  all  infomiation  is  specified  at  the  upstream  end.  The  example 
contains  a  transition  from  a  100-ft  (30.5-m)  width  to  a  90-ft  (27.4-m)  width. 
The  slope  is  0.01  except  in  the  lOO-ft  (30.5-m)  length  of  transition,  where  it  is 
0.02.  The  transition  is  to  have  a  bed  roughness  higher  than  the  rest  of  the 
model.  The  model  includes  a  50-ft  (15.2-m)-long  by  18-ft  (5.5-m)-wide  bridge 
pier.  Note  that  the  resolution  is  concentrated  around  and  downstream  of  any 
abrupt  changes  (Figure  C2).  Since  this  flow  is  supercritical,  all  information 
including  shocks  and  rough  conditions  is  swept  downstream.  Therefore,  this  is 
where  the  resolution  is  needed. 


Figure  C1.  Complete  example  grid  and  geometry 


Figure  C2.  Detail  of  example  grid  showing  transition  and  bridge  pier 
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The  complete  geometry  file  includes  338  nodes  and  302  elements.  The 
nodes  1,  2,  and  3  are  the  upstream  boundary  nodes  (designated  by  BI  in  the 
geometry  file),  and  nodes  333-338  are  the  downstream  boundary  nodes 
(designated  by  BO  in  the  geometry  file).  An  abbreviated  version  of  the 
geometry  file  is  shown  in  Figure  C3. 


T1 

example  geometry 

file. 

contains 
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Figure  C3.  Example  geometry  file  (example.geo) 


The  initial  hot  start  file  can  be  tricky.  This  case  involves  a  straight  reach 
so  it  is  fairly  easy  to  set  up  a  file  with  p,  q,  h  of  60., 0.,  3.  as  a  constant 
throughout  the  domain.  If  the  domain  is  curved  one  might  have  to  start  with 
pooled  flow  at  a  single  depth.  The  beginning  of  the  initial  hot  start  file  is 
shown  in  Figure  C4. 

The  initial  input  file  is  shown  in  Figure  C5. 

The  text  in  the  hydrodynamic  input  file  that  describes  each  numerical  input 
is  not  actually  read  by  the  program  so  it  can  be  left  out  or  changed.  The  ini¬ 
tial  time-step  was  chosen  to  produce  a  CFL  number  (Equation  15)  of  roughly  1 
(so  this  time-step  is  chosen  as  1.0  sec).  The  temporal  derivative 
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Figure  C4.  Example  hot  start  file  (example.hot). 
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Figure  C5.  Initial  hydrodynamic  input  file  (examplel.flo). 


was  chosen  to  be  first-order  since  the  interest  is  in  steady-state  results  only. 
Even  though  the  downstream  boundary  is  intended  to  be  supercritical  and  no 
boundary  condition  is  needed,  experience  has  shown  that  the  model  sometimes 
will  converge  to  a  different  solution  or  be  unstable  unless  a  very  good  starting 
guess  at  the  solution  was  made.  Therefore,  until  the  model  settles  down,  a 
tailwater  elevation  is  specified  corresponding  to  flow  that  is  slightly  subcri- 
tical.  Two  different  element  types  and  Manning's  roughness  coefficients  have 
been  designated.  The  convergent  section  is  steeper  and  rougher  in  this 
example.  The  model  will  run  20  time-stefs  resulting  in  a  simulation  time  of 
20  sec.  The  time-step  size  is  then  increased.  Later  the  downstream  tailwater 
condition  is  effectively  removed  by  specifying  an  elevation  that  is  below  the 
bed  elevation. 


The  next  input  file  is  run  to  raise  the  time-step  now  that  the  model  has 
settled  down  from  startup.  The  time-step  has  been  raised  to  3.0  sec.  The 
second  input  file  is  shown  in  Figure  C6. 


The  final  input  file  (Figure  Cl)  removes  the  tailwater  constraint.  This  run 
was  repeated  one  time. 


Figures  C8  and  C9  are  the  steady-state  results  of  the  model  run  for  water- 
surface  and  velocity  contours,  respectively.  These  were  generated  by  convert¬ 
ing  the  output  plot  file  to  a  binary  file  in  the  form  expected  for  a  solution  by 
FastTABS.  There  are  many  options  available  such  as  vectors  and  time-series 
plots.  The  results  show  that  the  bridge  pier  is  a  choke,  so  that  the  flow  is 
subcritical  in  front  of  the  pier.  This  is  apparent  since  the  shock  contour  is 
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Figure  C6.  Second  hydrodynamic  input  file  (example2.flo). 
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Figure  C7.  Final  hydrodynamic  input  file  (exampleS.flo). 


perpendicular  to  the  flow  instead  of  being  swept  back.  Note  also  that  the  flow 
returns  to  supercritical  flow  along  the  pier  sides  and  a  depression  wave  has  the 
sweptback  wake. 


Figure  C8.  Final  water-surface  contours,  ft 


Figure  C9.  Final  velocity  contours,  fps 
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Appendix  D 
Interpolation  Program 


The  interpolation  program  INTRPL8.EXE  allows  the  user  to  refine  a  grid 
for  which  answers  (p,  q,  and  h,  defined  in  Chapter  2)  have  already  been 
found.  Using  the  old  hot  start  information,  the  old  geometry,  and  the  new 
geometry,  the  user  can  generate  a  new  hot  start  file  for  the  new  geometry. 
This  new  hot  start  file  is  then  the  starting  point  for  a  new  run  of  HIVEL2D. 

The  program  also  outputs  an  interpolated  plot  output  file  that  can  be  con¬ 
verted  to  binary  and  viewed  in  FastTABS.  This  is  used  to  verify  the  inter¬ 
polated  values.  This  feature  can  also  be  used  to  interpolate  from  a  fine 
grid  to  a  coarse  grid,  to  better  view  vectors  on  a  very  fine  grid,  for  example. 
Users  should  note  that  the  program  cannot  interpolate  answers  for  any  nodes 
that  lie  outside  the  old  geometiy’s  boundary.  For  instance,  refining  a  mesh 
around  a  curve  will  cause  the  nodes  along  the  curve’s  outer  radius  to  lie 
outside  the  old  coarse  mesh’s  boundary.  The  program  will  report  the  nodes 
that  lie  outside  the  domain. 


Program  Description 

After  input  of  the  necessary  data,  the  first  routine  the  program  performs  is 
that  of  comparing  the  old  grid  with  the  new  grid.  Old  node  numbers  are 
typically  going  to  be  different  than  new  ones.  For  this  reason,  the  old 
answers  at  each  node  are  first  "mapped"  onto  the  corresponding  nodes  on  the 
new  grid  by  comparing  the  coordinates  of  the  nodes. 

Next  comes  the  task  of  looking  at  each  node,  seeing  if  it  is  a  new  node, 
and,  if  it  is  a  new  node,  finding  out  in  which  old  element  it  is  located.  There 
are  two  possible  locations  for  a  new  or  "target"  node:  in  a  quadrilateral  ele¬ 
ment  or  in  a  triangular  element.  A  node  on  an  element  side  is  interpolated  in 
the  element  in  which  it  is  first  found. 

Deciding  in  which  element  a  target  node  is  located  is  done  by  a  cross  prod¬ 
uct  scheme,  in  which  the  cross  product  is  taken  between  a  vector  lying  on  the 
side  of  an  element  and  a  vector  from  a  comer  node  on  that  side  to  the  target 
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point.  If  the  cross  products  for  all  sides  are  greater  than  zero,  then  the  node 
lies  within  that  element. 


Interpolation  Scheme  for  Quadrilateral  Elements 

For  convenience  in  interpolating,  a  local  coordinate  system  (^,11)  is  used 
where  -1  <  ^  <  1  and  -1  <  ti  <  1  (Figure  Dl)  to  derive  the  bilinear  interpola¬ 
tion  functions.  First  the  target  point  (Xp,yp)  is  written  in  (^,ri)  coordinates. 


Figure  D1 .  Local  coordinates  for  quadrilateral  and  triangular  elements 
This  is  done  using  Newton-Raphson  iteration  to  solve  the  following  equations: 


4 


4=1 


(Dl) 


and 


4 


y=T^iyi 

4=1 


(D2) 


for  the  local  coordinates  ^  and  q  where  the  bilinear  shape  functions  for  a 
quadrilateral  element  are: 


D2 
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<I>1  =  |(1-0(1-Tl) 

(i)2=l(i+oa-Ti) 

4)3=1(1  +  0(1+ti) 

4>4=Ta-oa+^) 


Interpolation  Scheme  for  Triangular  Elements 

Again,  the  local  coordinate  system  (^.Tj)  where  -1  <  ^  <  1  and  -1  <  t)  <1, 
(Figure  Dl)  is  used  to  derive  the  interpolation  functions.  Solving 


i=l 


y=E  Vi 

1=1 


simultaneously  for  \  and  T),  where: 


ioH) 


are  the  bilinear  shape  functions  for  a  triangular  element,  results  in  the 
following  equations  for  ^  and  rj: 


5  =  -  +  «2(-«3Pl  glP2) 

«i  «i(-Yi  Y2  Yj  -  Y4  -  Ys  Ye) 
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_ (-ttaPi  gip2) _  (D6) 

^  ("Yi  Y2  +  Yj  -  Y4  -  Ys  +  Ye) 

where 

«!  =  X,  -  Xj 

Ota  =  Xi  -  X3 

as  =  -  >2 

Pi  “  ~  “  •^s 

p2  =  -  )’2  -  ^3 

Yi  = 

Y2  =  ^3>1 

Y3  =  ^1^2 

Y4  =  ^3)'2 

Ys  = 

Y«  =  *2>’3 


Interpolation 

Once  the  target  node  has  been  located,  and  mapped  onto  old  grid  element 
local  coordinates  ^  and  T),  the  interpolated  values  of  p,  q,  and  h  are  found 
using  the  following  equations; 


'-'EM 


<=1 


(D7) 


1=1 


(D8) 
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where  n  is  the  number  of  nodes  on  the  old  grid  element  and  (j)  is  the 
appropriate  interpolation  function. 
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